options(future.globals.maxSize= 891289600)
### Make sure all the packages are installed
if (!requireNamespace("Seurat", quietly = TRUE))
install.packages("Seurat")
if (!requireNamespace("dplyr", quietly = TRUE))
install.packages("dplyr")
if (!requireNamespace("devtools", quietly = TRUE))
install.packages("devtools")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!requireNamespace("RcppML", quietly = TRUE)) {
# BiocManager::install("fgsea", update = FALSE)
# BiocManager::install("limma", update = FALSE)
# devtools::install_github("zdebruine/singlet", upgrade = FALSE)
devtools::install_github("zdebruine/RcppML")
devtools::install_github("zdebruine/RcppSparse")
}
if (!requireNamespace("sparseMatrixStats", quietly = TRUE)) {
install.packages("sparseMatrixStats")
}
if (!requireNamespace("inflection", quietly = TRUE)) {
install.packages("inflection")
}
if (!requireNamespace("org.Hs.eg.db", quietly = TRUE)) {
BiocManager::install("org.Hs.eg.db")
}
if (!requireNamespace("clusterProfiler", quietly = TRUE)) {
BiocManager::install("clusterProfiler")
}
### Load all the necessary libraries
library(Seurat)
library(dplyr)
library(RcppML)
library(RcppSparse)
library(ggplot2)
library(sparseMatrixStats)
library(inflection)
library(glue)
library(tidyr)
library(org.Hs.eg.db)
library(clusterProfiler)
set.seed(687)5 - Factor analysis
Introduction
In this notebooks we are going to carry out Factor Analysis analysis using RcppML. You can see the GitHub repository here. This implementation of NMF is extremely fast and enables us to use large dataset since it works well with sparse matrices. The author is currently implementing it within the Singlet package to work nicely with single cell and soon it will be available!
here and the differential communication analysis vignette here. The goal of the notebook is to identify which communication pathways are altered between flu and covid infection!
Some associated literature which is a must read are:
Why use NMF instead of PCA - This stats.StackExchange summarises quite well what the differences are between NMF and PCA. It uses the image below to visually represent the differences - see the free book An Introduction to Statistical Learning for more details!
Key Takeaways
NMF identifies sets of genes “metagenes” representing the main characteristics of the data.
Choosing the rank (K) of the matrix is an important step since it will determine how many “metagenes” are present in our dataset. We can use cross validation to find the best K.
Whe can visualize how important each genes is for each factor (aka - metagene) and how important is each factor for each cell. This way we can determine which metagenes are explaining each cell’s transcriptome.
When looking at metagenes learned across multiple patients and conditions it is imperative to check that the signal is not being provided just by one replicate of our experiment. Not checking this could lead to incorrect interpretation of the data!
Library
Load data
We’re going to be working with a dataset from the paper - Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19 Download data from cellxgene portal.
# Download the data in data/ directory
# download.file(
# url = "https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds",
# destfile = "../data/workshop-data.rds",
# method = "wget",
# extra = "-r -p --random-wait")
# We can also use the CLI with the wget command below
# wget https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds
se <- readRDS("../data/workshop-data.rds")Analysis
Convert ENSEMBL IDs to Gene Symbols
Right away we can see how ensembl ids are used in the rownames. Let’s transform them into their matched symbols to make them human-readable:
head(rownames(se))[1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457"
[5] "ENSG00000000460" "ENSG00000000938"
Convert to gene symbols
gene_df <- readr::read_csv(file = "../data/cov_flu_gene_names_table.csv")
symbol_id <- data.frame(index = rownames(se)) %>%
left_join(gene_df) %>%
pull(feature_name)
# re-create seurat object
mtx <- se@assays$RNA@data
rownames(mtx) <- symbol_id
se <- CreateSeuratObject(counts = mtx, meta.data = se@meta.data)# symbol_id <- mapIds(
# org.Hs.eg.db,
# keys = rownames(se),
# column = "SYMBOL",
# keytype = "ENSEMBL",
# multiVals = "first")
#
# # df <- data.frame(symbol = symbol_id, ensembl = names(symbol_id))
# all(rownames(se) == names(symbol_id))
#
# # re-create seurat object
# mtx <- se@assays$RNA@data
# rownames(mtx) <- symbol_id
#
# # Remove NAs
# sum(is.na(rownames(mtx)))
# dim(mtx)
# mtx <- mtx[!is.na(rownames(mtx)), ]
# dim(mtx)
#
# rownames(mtx) <- make.unique(rownames(mtx))
# se <- CreateSeuratObject(counts = mtx, meta.data = se@meta.data)
#
# sum(is.na(rownames(se)))Factor Analysis
For the purpose of this vignette we’re going to focus on CD8 and NK cells this object to pull out interferon signalling program:
se <- se[, se$Celltype %in% c("CD8, non-EM-like", "CD8, EM-like", "NK cell")]First we need to normalize the data
se <- NormalizeData(se, verbose = FALSE) %>%
FindVariableFeatures(nfeatures = 3000, verbose = FALSE) %>%
ScaleData(verbose = FALSE) %>%
RunPCA(verbose = FALSE)
VlnPlot(
se,
features = c("CD3D", "CD3E", "CD8B", "NKG7", "FCGR3A"),
group.by = "Celltype",
pt.size = 0,
split.by = "disease") + theme(legend.position = "bottom")Extract the normalized expression matrix and remove genes that are all 0s
library(Matrix)
mtx <- se@assays$RNA@layers$data
colnames(mtx) <- colnames(se)
rownames(mtx) <- rownames(se)
# Remove genes that are all 0s
mtx <- mtx[sparseMatrixStats::rowSums2(mtx) > 0, ]Cross-validation for rank determination
The first step is to determine the rank of our matrix - by this we mean, which is the optimal number of factors we need to decompose our matrix. Here we run cross-validation 3 times across a range of different ranks (k).
start <- Sys.time()
cv_data <- crossValidate(mtx, k = seq(1, 36, 5), reps = 2, verbose = TRUE, tol = 1e-04)
Sys.time() - startTime difference of 3.998327 mins
Now let’s visualize the crossvalidation results
ggplot(cv_data, aes(x = k, y = value, color = rep, group = rep)) +
geom_point() +
geom_line() +
theme_minimal()We can see how it starts to plateau at ~k=20 so we will use that in this script.
Run NMF
k: Number of factors we want to identify
tol: Stands for tolerance and is how small we want the error to be, the smaller the number the better the decomposition but the longer its gonna take to converge. Values in the 10^-5 and 10^-6 return very good results.
L1: Introduces sparsity into the factors and loadings with the aim of removing noisy genes from factors and noisy cells contributing to factors. L1 normalization uses Lasso regularization to increase sparsity and remove the unimportant features.
set.seed(7)
nmf_ls <- RcppML::nmf(
data = mtx,
k = 20,
tol = 1e-06,
L1 = c(0.05, 0.1), # l1 for c(w, h)
verbose = TRUE)Explore the NMF results
W - contains which genes are relevant for each factor. H - contains how important each factor for each cell.
w <- nmf_ls@w
w[1:5, 1:5] nmf1 nmf2 nmf3 nmf4 nmf5
TSPAN6 0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
DPM1 0 9.446519e-05 4.755987e-05 1.160628e-04 7.280245e-05
SCYL3 0 0.000000e+00 0.000000e+00 5.235962e-07 0.000000e+00
C1orf112 0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
FGR 0 5.544984e-05 4.062468e-04 5.537400e-05 2.077896e-04
h <- nmf_ls@h
h[1:5, 1:5] AAACCCAAGAATGTTG-12 AAACCCAAGCTACGTT-6 AAACCCAAGGCCTGCT-12
nmf1 5.462850e-05 5.187584e-05 3.464711e-05
nmf2 7.132754e-05 1.925005e-04 0.000000e+00
nmf3 1.210247e-04 2.170797e-05 0.000000e+00
nmf4 1.459895e-04 0.000000e+00 0.000000e+00
nmf5 5.263035e-05 0.000000e+00 1.276532e-04
AAACCCAAGTGTAGAT-5 AAACCCACACAATGCT-5
nmf1 6.166142e-05 5.885350e-05
nmf2 9.268590e-05 1.126141e-04
nmf3 0.000000e+00 0.000000e+00
nmf4 4.154653e-06 2.351866e-05
nmf5 7.246126e-05 1.722137e-05
Add NMF to Seurat object
se[["FA"]] <- Seurat::CreateDimReducObject(
embeddings = t(nmf_ls$h),
loadings = nmf_ls$w,
assay = "RNA",
key = "Factor_")
se[["FA"]]A dimensional reduction object with key Factor_
Number of dimensions: 20
Number of cells: 19262
Projected dimensional reduction calculated: FALSE
Jackstraw run: FALSE
Computed using assay: RNA
Examine and visualize NMF loadings a few different ways
print(se[["FA"]], dims = 1:20, nfeatures = 5)Factor_ 1
Positive: MALAT1, MT-CO3, MT-CO2, RPL10, EEF1A1
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 2
Positive: GNLY, CD52, GZMH, FGFBP2, NKG7
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 3
Positive: GNLY, GZMB, NKG7, GZMH, IL32
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 4
Positive: GNLY, FGFBP2, NKG7, TYROBP, CCL5
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 5
Positive: NFKBIA, GNLY, CD69, CCL4, IER2
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 6
Positive: JUN, JUNB, CXCR4, DUSP1, TSC22D3
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 7
Positive: S100A9, S100A8, PRF1, NKG7, GZMA
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 8
Positive: CMC1, CD3D, NKG7, HLA-DRB1, CD74
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 9
Positive: ACTB, LGALS1, CALR, HSP90AB1, PFN1
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 10
Positive: MALAT1, MT-CO2, MT-CO3, MT-CO1, MT-ATP6
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 11
Positive: CREM, METRNL, CXCR4, DUSP2, PIK3R1
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 12
Positive: FCER1G, SPON2, MYOM2, FGFBP2, IGFBP7
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 13
Positive: CREM, ZNF331, AREG, CEMIP2, DUSP2
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 14
Positive: CD74, ACTG1, HLA-DRA, COTL1, LIMD2
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 15
Positive: GNLY, CTSW, SELL, GZMK, CD7
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 16
Positive: ISG15, MTRNR2L12, XAF1, JUN, IFI6
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 17
Positive: LTB, IL7R, CD52, VIM, IL32
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 18
Positive: RPS26, METRNL, HLA-DQB1, HLA-DRB1, ID2
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 19
Positive: GNLY, GZMB, PRF1, GZMH, FLNA
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 20
Positive: KLRD1, TXNIP, ZFP36L2, KLF2, RPS4Y1
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Violin plots for the loadings
VlnPlot(
se,
features = glue::glue("Factor_{1:ncol(nmf_ls$w)}"),
ncol = 4,
group.by = "Celltype")We can also visualize this by cell type, condition and patient!
disease_pal <- c("#41AE76", "#225EA8", "#E31A1C")
names(disease_pal) <- c("normal", "influenza", "COVID-19")
# flu <- RColorBrewer::brewer.pal(12, name = "YlGnBu")
# normal <- RColorBrewer::brewer.pal(12, name = "BuGn")
# covid <- RColorBrewer::brewer.pal(12, name = "YlOrRd")
donor_pal <- c(
"#66C2A4", "#41AE76", "#238B45", "#006D2C",
"#41B6C4", "#1D91C0", "#225EA8", "#253494", "#081D58",
"#FFFFCC", "#FFEDA0", "#FED976", "#FEB24C", "#FD8D3C",
"#FC4E2A", "#E31A1C", "#BD0026", "#800026")
names(donor_pal) <- c(
"Normal 1", "Normal 2", "Normal 3", "Normal 4",
"Flu 1", "Flu 2", "Flu 3", "Flu 4", "Flu 5",
"nCoV 1", "nCoV 2", "nCoV 3_4", "nCoV 5",
"nCoV 6", "nCoV 7_8", "nCoV 9_10", "nCoV 11"
)
# Preprocess dataset
dd <- bind_cols(se@meta.data, se@reductions$FA@cell.embeddings) %>%
tidyr::pivot_longer(
cols = glue::glue("Factor_{1:ncol(nmf_ls$w)}"),
names_to = "k",
values_to = "loading") %>%
group_by(Celltype, disease, donor_id, k) %>%
summarise(median_loading = median(loading))
lapply(glue::glue("Factor_{1:ncol(nmf_ls$w)}"), function(i) {
dd %>%
filter(k == i) %>%
ggplot(aes(x = disease, y = median_loading, fill = disease)) +
geom_boxplot(outlier.shape = NA, alpha = 0.7) +
geom_jitter(aes(color = donor_id), height = 0) +
facet_wrap(~ Celltype) +
scale_fill_manual(values = disease_pal) +
scale_color_manual(values = donor_pal) +
scale_x_discrete(limits = names(disease_pal)) +
theme_classic() +
labs(title = glue::glue("Factors by condition in {i}"))
})[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
[[9]]
[[10]]
[[11]]
[[12]]
[[13]]
[[14]]
[[15]]
[[16]]
[[17]]
[[18]]
[[19]]
[[20]]
Factors 4 & 6 seem interesting, lets explore them a bit more in depth!
w[1:5, c("nmf5", "nmf7")] nmf5 nmf7
TSPAN6 0.000000e+00 0.0000000000
DPM1 7.280245e-05 0.0000000000
SCYL3 0.000000e+00 0.0000000000
C1orf112 0.000000e+00 0.0000000000
FGR 2.077896e-04 0.0005231307
# Extract top genes from factors 4 & 6 using Unit Invariant Knee
FA_genes <- lapply(c("nmf5", "nmf7"), function(i) {
y_vec <- sort(w[, i], decreasing = TRUE)
# Define inflecion point
n <- uik(x = seq_len(length(y_vec)), y = y_vec)
# Define top 10 genes
df <- data.frame(gene = names(y_vec), value = y_vec, rank = 1:length(y_vec)) %>%
mutate(lab = if_else(rank < 15, gene, NA_character_))
print(ggplot(df, aes(x = rank, y = y_vec)) +
geom_point() +
geom_vline(xintercept = n) +
ggrepel::geom_text_repel(aes(label = lab)) +
labs(title = glue("Factor-{i}")) +
theme_classic())
y_vec[1:n]
})names(FA_genes) <- c("nmf5", "nmf7")
dd <- lapply(names(FA_genes), function(i){
data.frame(
value = FA_genes[[i]],
gene = names(FA_genes[[i]]),
factor = i)
}) %>%
bind_rows() %>%
pivot_wider(names_from = factor, values_from = value, values_fill = 0)
DT::datatable(dd)By looking at these genes we can get a sense of which are the major processes captures by each factor.
NMF-5 seems to be capturing cytotoxic and activation signals since it contains genes encoding for chemokines (CCL4, CCL3), granzymes (GZMB, GNLY, NKG7) and other immune related processes (NFKBIA, CD69…).
NMF-7 contains a wide array of myeloid cell markers, such as S100A8, S100A9, LYZ and genes encoding for MHC-II (CD74, HLA-DRB1, HLA-DRA) as well as cytotoxic genes (GZMA, GZMH, NKG7…). The most straight forward explanation is that these could be doublets or have a high amount of ambient RNA in the soup.
We can take a look at how these genes are expressed across the patients since like Flu-1 patient has an etremely high score and we want to make sure it is not driving this by itself.
# We'll take a look at the top 25 genes
names(head(FA_genes[["nmf7"]], 25)) [1] "S100A9" "S100A8" "PRF1" "NKG7" "GZMA" "IL32" "CD74"
[8] "GZMB" "LYZ" "GZMH" "HLA-B" "B2M" "HBB" "TXNIP"
[15] "S100A11" "MT-CO1" "PSME2" "PSMB9" "PLAAT4" "RPS4Y1" "CST7"
[22] "TRAC" "FGFBP2" "ACTB" "HLA-A"
# Scale their expression
se <- ScaleData(se, features = names(head(FA_genes[["nmf7"]], 25)))
# Check scale.data
mtx <- se@assays$RNA@layers$scale.data
rownames(mtx) <- names(head(FA_genes[["nmf7"]], 25))
colnames(mtx) <- colnames(se)
# Choose color palette
palette_length <- 100
my_color <- colorRampPalette(c("Darkblue", "white","red"))(palette_length)
# set color breaks
my_breaks <- c(seq(quantile(mtx, .01), 0, length.out=ceiling(palette_length/2) + 1),
seq(0.05, quantile(mtx, .99), length.out=floor(palette_length/2)))
# Note that we are setting the min and max of the scale to .01 and .99 so we exclude extreme values from dampening our signal
# We can see how the max value of the matrix is +10 but q99 is 2
Hmisc::describe(as.vector(mtx))as.vector(mtx)
n missing distinct Info Mean Gmd .05
481550 0 287115 1 -2.226e-06 1.1 -1.59884
.10 .25 .50 .75 .90 .95
-1.31044 -0.59795 0.02677 0.69204 1.16664 1.48214
lowest : -12.5361 -8.98419 -8.74596 -8.11004 -7.98887
highest: 9.87275 9.88233 9.96698 9.97462 10
pheatmap::pheatmap(
mtx,
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = TRUE,
show_colnames = FALSE,
treeheight_col = 0,
annotation_col = se@meta.data[, c("disease", "donor_id", "Celltype")],
annotation_colors = list("disease" = disease_pal, "donor_id" = donor_pal),
color = my_color,
breaks = my_breaks)Carry out GSEA
# read GSEA markers
gsea_ls <- lapply(FA_genes, function(i) {
# http://yulab-smu.top/clusterProfiler-book/chapter5.html#go-gene-set-enrichment-analysis
gsea_results <- clusterProfiler::gseGO(
geneList = i,
ont = "BP",
OrgDb = org.Hs.eg.db,
keyType = "SYMBOL",
minGSSize = 10,
maxGSSize = 300,
pvalueCutoff = 0.1,
pAdjustMethod = "BH",
seed = TRUE)
gsea_results
})
names(gsea_ls) <- names(FA_genes)Visualize top enriched gene sets per cluster
lapply(names(gsea_ls), function(i) {
# Extract gsea
gsea <- gsea_ls[[i]]
gsea <- clusterProfiler::simplify(gsea, cutoff = 0.7)
gsea@result <- gsea@result %>%
dplyr::arrange(dplyr::desc(NES)) %>%
dplyr::filter(p.adjust < 0.1)
tmp_plt <- gsea@result %>%
dplyr::top_n(n = 20, wt = enrichmentScore) %>%
ggplot2::ggplot(.,
ggplot2::aes(
x = NES,
y = forcats::fct_reorder(Description, NES),
size = setSize,
color = p.adjust)) +
ggplot2::geom_point() +
# ggplot2::geom_vline(
# xintercept = 0,
# linetype = "dashed",
# color = "red") +
scale_color_viridis_c(option = "plasma") +
ggplot2::theme_minimal() +
ggplot2::labs(title = i)
}) %>% patchwork::wrap_plots(ncol = 2)And lastly we will save GSEA to a spreadsheet so we can take a look:
gsea_xlsx <- lapply(names(gsea_ls), function(i) {
print(i)
# Extract gsea
gsea <- gsea_ls[[i]]
gsea <- clusterProfiler::simplify(gsea, cutoff = 0.7)
gsea@result <- gsea@result %>%
dplyr::arrange(dplyr::desc(NES)) %>%
dplyr::filter(p.adjust < 0.1)
})[1] "nmf5"
[1] "nmf7"
names(gsea_xlsx) <- names(gsea_ls)
openxlsx::write.xlsx(gsea_xlsx, file = "../data/GSEA.xlsx")Session Info
sessionInfo()R version 4.3.1 (2023-06-16)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] Matrix_1.6-5 clusterProfiler_4.8.2 org.Hs.eg.db_3.17.0
[4] AnnotationDbi_1.62.2 IRanges_2.34.1 S4Vectors_0.38.1
[7] Biobase_2.60.0 BiocGenerics_0.46.0 tidyr_1.3.1
[10] glue_1.7.0 inflection_1.3.6 sparseMatrixStats_1.12.2
[13] MatrixGenerics_1.12.3 matrixStats_1.2.0 ggplot2_3.5.0
[16] RcppSparse_1.0 RcppML_0.5.6 dplyr_1.1.4
[19] Seurat_5.0.1 SeuratObject_5.0.1 sp_2.1-3
loaded via a namespace (and not attached):
[1] fs_1.6.3 spatstat.sparse_3.0-3 bitops_1.0-7
[4] enrichplot_1.20.0 devtools_2.4.5 HDO.db_0.99.1
[7] httr_1.4.7 RColorBrewer_1.1-3 backports_1.4.1
[10] profvis_0.3.8 tools_4.3.1 sctransform_0.4.1
[13] DT_0.29 utf8_1.2.4 R6_2.5.1
[16] lazyeval_0.2.2 uwot_0.1.16 urlchecker_1.0.1
[19] withr_3.0.0 prettyunits_1.1.1 gridExtra_2.3
[22] progressr_0.14.0 cli_3.6.2 spatstat.explore_3.2-6
[25] fastDummies_1.7.3 scatterpie_0.2.1 sass_0.4.8
[28] labeling_0.4.3 spatstat.data_3.0-4 readr_2.1.4
[31] ggridges_0.5.6 pbapply_1.7-2 yulab.utils_0.1.4
[34] foreign_0.8-84 gson_0.1.0 DOSE_3.26.2
[37] parallelly_1.37.0 sessioninfo_1.2.2 rstudioapi_0.15.0
[40] RSQLite_2.3.1 generics_0.1.3 gridGraphics_0.5-1
[43] crosstalk_1.2.1 ica_1.0-3 spatstat.random_3.2-2
[46] vroom_1.6.3 zip_2.3.0 GO.db_3.17.0
[49] ggbeeswarm_0.7.2 fansi_1.0.6 abind_1.4-5
[52] lifecycle_1.0.4 yaml_2.3.8 qvalue_2.32.0
[55] Rtsne_0.17 grid_4.3.1 blob_1.2.4
[58] promises_1.2.1 crayon_1.5.2 miniUI_0.1.1.1
[61] lattice_0.21-8 cowplot_1.1.3 KEGGREST_1.40.0
[64] pillar_1.9.0 knitr_1.45 fgsea_1.26.0
[67] future.apply_1.11.1 codetools_0.2-19 fastmatch_1.1-4
[70] leiden_0.4.3.1 downloader_0.4 ggfun_0.1.4
[73] data.table_1.15.0 remotes_2.4.2.1 vctrs_0.6.5
[76] png_0.1-8 treeio_1.24.3 spam_2.10-0
[79] gtable_0.3.4 cachem_1.0.8 openxlsx_4.2.5.2
[82] xfun_0.42 mime_0.12 tidygraph_1.2.3
[85] survival_3.5-7 pheatmap_1.0.12 ellipsis_0.3.2
[88] fitdistrplus_1.1-11 ROCR_1.0-11 nlme_3.1-163
[91] ggtree_3.8.2 usethis_2.2.2 bit64_4.0.5
[94] RcppAnnoy_0.0.22 GenomeInfoDb_1.36.3 bslib_0.6.1
[97] irlba_2.3.5.1 rpart_4.1.19 vipor_0.4.5
[100] KernSmooth_2.23-22 Hmisc_5.1-0 colorspace_2.1-0
[103] DBI_1.1.3 nnet_7.3-19 ggrastr_1.0.2
[106] tidyselect_1.2.0 processx_3.8.2 bit_4.0.5
[109] compiler_4.3.1 htmlTable_2.4.1 plotly_4.10.4
[112] shadowtext_0.1.3 checkmate_2.2.0 scales_1.3.0
[115] lmtest_0.9-40 callr_3.7.3 stringr_1.5.1
[118] digest_0.6.34 goftest_1.2-3 spatstat.utils_3.0-4
[121] rmarkdown_2.25 XVector_0.40.0 base64enc_0.1-3
[124] htmltools_0.5.7 pkgconfig_2.0.3 fastmap_1.1.1
[127] rlang_1.1.3 htmlwidgets_1.6.4 shiny_1.8.0
[130] jquerylib_0.1.4 farver_2.1.1 zoo_1.8-12
[133] jsonlite_1.8.8 BiocParallel_1.34.2 GOSemSim_2.26.1
[136] RCurl_1.98-1.12 magrittr_2.0.3 Formula_1.2-5
[139] GenomeInfoDbData_1.2.10 ggplotify_0.1.2 dotCall64_1.1-1
[142] patchwork_1.2.0 munsell_0.5.0 Rcpp_1.0.12
[145] ape_5.7-1 viridis_0.6.4 reticulate_1.35.0.9000
[148] stringi_1.8.3 ggraph_2.1.0 zlibbioc_1.46.0
[151] MASS_7.3-60 plyr_1.8.9 pkgbuild_1.4.2
[154] parallel_4.3.1 listenv_0.9.1 ggrepel_0.9.5
[157] forcats_1.0.0 deldir_2.0-2 Biostrings_2.68.1
[160] graphlayouts_1.0.0 splines_4.3.1 tensor_1.5
[163] hms_1.1.3 ps_1.7.5 igraph_2.0.2
[166] spatstat.geom_3.2-8 RcppHNSW_0.6.0 reshape2_1.4.4
[169] pkgload_1.3.2.1 evaluate_0.23 BiocManager_1.30.22
[172] tzdb_0.4.0 tweenr_2.0.2 httpuv_1.6.14
[175] RANN_2.6.1 purrr_1.0.2 polyclip_1.10-6
[178] future_1.33.1 scattermore_1.2 ggforce_0.4.1
[181] xtable_1.8-4 RSpectra_0.16-1 tidytree_0.4.6
[184] later_1.3.2 viridisLite_0.4.2 tibble_3.2.1
[187] aplot_0.2.2 memoise_2.0.1 beeswarm_0.4.0
[190] cluster_2.1.4 globals_0.16.2